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ABSTRACT 

Baryon acoustic oscillations are an excellent technique to constrain the properties of dark energy 
in the Universe. In order to accurately characterize the dark energy equation of state, we must 
understand the effects of both the nonlinearities and redshift space distortions on the location and 
shape of the acoustic peak. In this paper, we consider these effects using the Zel'dovich approximation 
and a novel approach to 2nd order perturbation theory. The second order term of the Zel'dovich power 
spectrum is built from convolutions of the linear power spectrum with polynomial kernels in Fourier 
space, suggesting that the corresponding term of the the Zel'dovich correlation function can be written 
as a sum of quadratic products of a broader class of correlation functions, expressed through simple 
spherical Bessel transforms of the linear power spectrum. We show how to systematically perform 
such a computation. We explicitly prove that our result is the Fourier transform of the Zel'dovich 
power spectrum, and compare our expressions to numerical simulations. Finally, we highlight the 
advantages of writing the nonlinear expansion in configuration space, as this calculation is easily 
extended to redshift space, and the higher order terms are mathematically simpler than their Fourier 
counterparts. 

Subject headings: cosmology: theory - large-scale structure of Universe 



1. INTRODUCTION 

Baryon acoustic oscillations (BAO) arise from photon pressure opposing the g ravitational collapse of baryonic density 
perturbations in the baryon-photon plasma be fore the epoch of recombination ( Peebles fc Yu|1970 Silk|1968 Sunyaev 



1970 Bond & Efstathiou 1984). We se e the signature of this ettect m both the^ tem perature power 
le cosmic microwave background (CMB) ( de Bernardis et aL||2000 'Hanany et al.|[2000 |, and the galaxy 



fc Zel'dovich 

spectrum of t 

power spectrum. The BAO feature was first detected as a peak in the correlation function of the Luminous Red 
Galaxies (LRGs) in the Sloan D igital Sky Survey (SPSS) by Eisenstein et al. ( 20 05j), and was subsequ ently measured 
in several other galaxy surveys ( Cole et ar]|2005 Padmanabhan et al.|[2b07j Percival et al. [2007 2010). 

The BAO signal is useful because at low redshift, the scale of interest is still largely in the linear regime. The location 
of the peak thus acts as a "standard ruler" with which the expansion history of t he universe can b e measured as a 
function of time. This allows us to constrain the equation of state of Dark Energy ( Seo fc Eisenstein|[200 3 ) . 

In order to accurately constrain the properties of dark energy using BAO, we must understand the effects of non- 
linearities and redshift space distortions on the acoustic peak. Linear theory predicts that the amplitude of the peak 
will scale as the square of the linear growth function, D(t). However, it is important to understand to what extent 
the nonlinearities contribute to the location and shape of the peak at low redshift. It is also well kn own that there are 
first order effects on the correlation function due to the transformation from real to redshift space ( |Kaiser|[T987| [Tia^ 
et al.|[20TT Blake et al.||2011 Seo et al~]|2012 1 . Understanding how this transformation affects higher order terrns will 
be increasingly important as observations become more precise. 
Much work has been done to characterize the nonlinearities using cosmological perturbation theory dJain & 




many advantages over the power spectrum when analyzing the HAU sig: 
a method of calculating the first nonlinear term in the real-space correlation function using the Zel'dovich approxima- 
tion. In Section [2] we show the derivation of the nonlinear term from the Zel'dovich approximation. In Section [3] we 
show how this is the exact Fourier transform of the nonlinear Zel'dovich power spectrum. We dir ectly compare our 
calcul ation to the expressio n for the nonline ar Zel'dovich power spectrum, originally calculated in Grinstein & Wise 
(1987), and more recently in Valageas (2011). We show in Section B] that our configuration-space term is in agreement 
with numerical simulations of the Zel dovich approximation. We conclude in Section [5] 



2. COMPUTATION OF THE NONLINEAR ZEL'DOVICH CORRELATION FUNCTION IN REAL SPACE 

The Zel'dovich approximation maps particles' initial Lagrangian coordinate, q, to their co-moving Eulerian coordi- 
nate, x, through the gradient of the linear potential, (^(q), and the growth function D{t) ( Zerdovicli||1970 Shandarin 



fc Zeldovich||1989 ). The linear potential is related to the linear density field through the Poisson equation: 

x(q,t) =q-i?(t)V,0(q) , 
i?(i)V20(q)=<5i(q,i) . 



(1) 



In the rest of this paper, we will shorten D{t) to D for convenience. The real-space density is then related to the 
Jacobian of the transformation from Lagrangian to Eulerian coordinates. 



dxi 



1 



J(q,t) 



= 1 + 5(q(x)) , 



(2) 



where x and q are related by the Zel'dovich formula ([T]). Here, S (without subscript L) is the weakly nonlinear over- 
density. Equation (pi) for the Euleria n density is on l y stric tly valid before shell-crossing, where the mapping from q 
to X is one-to-one. As is discussed in Kofman et al.| ( |1994| , when there is multi-streaming, multiple values of q = q, 
map to the same point x, thus the Eulerian density at a given point is a sum over all of the streams at that point. 



P(x) 



E 

1=1 



(3) 



Here, p(x) is the full Eulerian density, and g{q_i) are the "single-stream" densities. However, shell-crossing does not 
contribute greatly at the large scales we are concerned with, and so we assume the "single-stream" density is the full 

density. We will show that this assumption leads to a result that agrees with the usual perturbative results. 

The Jacobian can be written in terms of invariants, /i, I2, and I3, of the symmetric matrix dij{q) ( Zel ' dovich| 1970 ). 



rfjj(q) 



920(q) 



(4) 
(5) 



dqidqj 

J(q, t) = 1 - Dh (q) + D^h (q) - D^h (q) ■ 
The invariants can be written in terms of the eigenvalues of dij(q): Ai, A2, and A3. 

/i(q) = Ai A2 -I- A3 , 

^2(q) = A1A2 A1A3 + A2A3 , (6) 
4(q) = A1A2A3 . 

Therefore, the Eulerian overdensity can be expressed by the Taylor expansion of the inverse Jacobian, to any order: 

6iq,t) = Dhiq) + (hiqf - h{q)) + {h{qf - 2/i(q)/2(q) + /3(q)) + ... (7) 

Note that the first term in this expansion corresponds to linear theory: 6l = DIi{q). 

However, as our goal in this paper is to calculate the density correlation function at a fixed Eulerian distance, we 
need to be careful about the details of the difference between the Eulerian and Lagrangian coordinates. Instead of 
the usual forward relation, writing x in terms of q, we will need the reverse and express q with x, using the Taylor 
expansion around x in a recursive fashion. 



q(x) =x + i?V5(/)(q(x)) 



(8) 



To zeroth order, Vg(/)(q(x)) = Vg(/i(x), but when considering the higher order terms in the density, it is necessary to 
express \/q(j){q{x)) to linear order in D. 

dH^iq)) 



dqi 



dm 
dq^ 



dcj^jq) 
dqidqj dq^ 



(9) 



q=x 



Because the correlation function, ^(r), is a function of Eulerian distance, r = Xi — X2, we must express the over- 
density as a function of the Eulerian coordinate, x. By expanding the function 6{q) about the point q = x, we arrive 
at an expression for (5(x) that is a power series in D. 



^(x,O='5(x + I?V,0(q(x))) . 



A Taylor expansion of the right hand side of Equation (|10|) gives 
,5(x,t)= (d{q,t) + DY, 



90(q(x))95(q,t) 



dqi 



dqi 



d^S{q,t) dcf^jq) 50(q) 
dqidqj dqi dqj 



(10) 



(11) 



With this, we can write the Eulerian density in terms of the Lagrangian density to third order in D: 

,^ f,f„_,, , ^s:^ d(j){ci) dS{q,t) 2X^ 9^<l>{<l)d(l){q)dS{q,t) 
Si^,t)^[Si<,,t)+D^—^—+D La^^^^ 

2 ^ dqidqj dqi dq. 



(12) 

q=x 

We use Equation ([?]) to express 5 in terms of the hnear quantities 6l and 0, where 5l is already first order in D. 
The real-space correlation function in co-moving Eulerian coordinates is: 

e(xi-X2,t) = (5(xi,t)(5(x2,t)) . (13) 

This can be written in powers of D using the above expansion of the Eulerian over-density. Because the over-density 
field is assumed to be a zero-mean Gaussian random field, the odd moments vanish. The correlation function to second 
order is then: 

We define the functions: 

C(0 = / PL(A:)j„(fcr)fc"+2dfc , (14) 

where j„ is the spherical Bessel function of order n and PL{k) is the linear power spectrum. Using this definition, the 
linear term is: 

e^'Hr)-eSw, (15) 

the spherically symmetric Fourier transform of the linear power spectrum. 

(r) is the expectation of a sum of products of four terms expressed with the linear quantities. Since these are 
all Gaussian, the only irreducible terms are second order. Thus we only need to calculate expectations of the type 
(a(qi)5(q2)). Mathematica was used to express the various derivatives of ^(q) in terms of spherical harmonics, and 
to calculate the expectation values between them. We illustrate this process in the following example. 



Oq. - y (2.)3fc.^-«^'''''--VTy(2^^^^«^^'''^' (16) 
^''^(q) _ f d^l^^^^^yy..^ fd^ f^,^^ ^ ^^^^ A.^V^.^k.q 



^ (Yo'ie,^) + ^Y°{e,^)j SUk)e^^-^ . (17) 



(27r)3/c2 ' 3 7 (2 

/ a0(qi) 920(q2) \ _ ini f [ {5L{\^i)h(}^2)) 



\ dq, dql I 33/2(27r)6 

Xy°((?i,^l) (^y„°(02,¥'2) + ^i^2°(^2,'/'2)) . (18) 

(<5L(ki),5L(k2)) - (2^)3,5^(ki + k2)PL(fcl) . (19) 

So we must calculate the integrals 

V,"(6',95)y,r"'(6i,v3)e**'-(ii-'i^)df7fc . (20) 



We use the plane- wave (Rayleigh) expansion for e^^''^^'^ i^) to compute these integrals in terms of spherical Bessel 
functions, jnikr), where r = qi — q2. In general, the integral is zero for m ^ m' . 



\ dq, dql I 5(2^2) 



PL{k) {3ji{kr)Vi{fi) - 2j3{kr)r3{fi))kdk , 







= -l'Pi{f^)^^\r) + lv3{fi)^3\r) . (21) 

where is the cosine of the angle between f and z, and Vi{fJ.) are the Legendre polynomials. In the expansion of the 
over-density. Equation (12 1, each term is evaluated at q = x, so the expectation values in the correlation function 
become functions of theTTulerian distance, r = |xi — X2|. 



4 



In the full calculation, there are many seemingly serendipitous cancellations of terms, which in the end give us a 
relatively simple expression. As we found out, any small mistake in the Mathematica script can lead to a far more 
complicated expression with the wrong asymptotics, so it is important to understand the expected behavior of the full 
correlation function to check the result. Firstly, the expectation value of the over-density should vanish to second order. 
Secondly, the correlation function should be isotropic in real space. Finally, we consider the asymptotic behavior of 
the expression - at zero lag, it must be finite, as this is the variance of the over-density field, and on large scales, the 
expression must tend to zero. Requiring this behavior rules out many erroneous results, and was a useful debugging 
tool. 

The full expression for the nonlinear Zel'dovich correlation function in Eulerian coordinates is: 

- ^Ca-'MeKO + lio'ir)eoir) + ^e2"'(0el(0 • (22) 

?(^)(0) = geoW = |ao^ (23) 

where ctq is the variance of the linear over-density field. As expected, the nonlinear term is finite at r = 0. As r — > oo, 
each set of terms within the expression tends to zero, due to the behavior of the spherical Bessel functions, so the full 
expression tends to zero as expected. Note that this expression contains products of the functions Ci?{r). In the next 
section, we show that this is the exact Fourier transform of the nonlinear Zel'dovich power spectrum. 

3. NONLINEAR POWER SPECTRUM 
The n ormalization we use for the power spectrum requires that we rewrite the nonlinear Zel'dovich power spectrum 



At zero lag, this takes the value: 



given in Valageas (20111 as: 



= -e<jlPL(k) 
1 



(24) 



fhfhx 1, , (k-ki)2(k-k2)^ ^ ^ 

(2^)3 '^i^(k-ki-k2) PL{kl)PL{k2) , 



PL{w)dw . 



We group Equation (22) into two terms, in order to clearly relate it to Piz(k) and P22(k): 

?^'nO = C(i3)('^)+C(22)(r) , (25) 

where C(i3)(f') is defined as the first term in Equation ( |22[ ), and 'C(22)('") contains the remaining terms. The term 
C(i3)('') straightforwardly Fourier transforms to Pi^{k): 



3(27r) 



lCo"'(0) 



OO nOO 



Pl {k')k"^io {k'r)jo {krydk'dr 



(26) 



^0-^(0)/ p^{k')k"' — Soik-k')dk 



Stt 



= -:;^om'PLik) 



This is equal to Pi3{k) from Equation (l24 



Pi,{k) = -k^alPLik) = --i-^{Q)k^PL{k) . 
The Fourier transform of C{22)(^) contains integrals of three spherical Bessel functions. In general: 

-^{CKOCrMI-^y^ PL(fci)PL(fc2)A;r^+'fcr+'dfcidfc2 ^ 3n{kir)j„{k2r)jo{krydr 
We require the expressions for integrals of the type 

o 

jo{kr)jn{kir)jn{k2r)r'^dr, 



(27) 



(28) 



(29) 



which can be found in Mehrem et al. ( 1991 ). 
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For example, to compute -7^{C2(^)^}- 

r Mk^r)Mk,r)jo{krydr = "-^^^^ m + ^1 - k'f ~ Aklkl) , (30) 

{1 \ki " k2\ < k < ki + k2 
\ \ki — k2\ ~ k or ki + k2 = k 
otherwise 

This means that the integral is zero unless k, ki, and k2 form a triangle. In the fci k2 plane, we consider the region 
where this triangular condition is satisfied. On the boundary of that region, the Bessel integral is discontinuous. 
However, since the boundary has measure zero, it does not contribute to the integral. Thus, we consider only the 
inside of the region of the plane where k, ki, and k2 form a triangle, i.e. where (3 — 1. 
In the integral, we make the substitution /c2 = |k — ki|, and integrate over the angle 6 between k and ki. 



m°Mn = ^ /d'fcli^L(fcl)PL(|k-ki| 



3((k - ki)^ + fcf - k'^ f - 4fcf (k - ki)2 

' fc?|k-ki|2 
^^^^^,jj<i'k,6i^k2PL{k^)PL{k2)5D{\^-\^l-\^2 r ' g^2 ' ' ' (31) 

This illustrates how the products of two in configuration space correspond to a convolution in Fourier space. 

The Fourier transform of all the terms in '^(22) ('') is: 

HMr)}- J f '^^^Snik - ki - k2) • "^^ff; PUki)PL{k2) . (32) 

which is exactly the P22{k) term from Equation (24 1. Thus, the entire expression transforms to P^'^\k). 

This agreement validates our configuration-space approach to perturbation theory with the Zel'dovich approximation. 
It is also remarkable how the triangular closure of the three wave vectors is connected to the integral over three spherical 
Bessel functions. 

4. COMPARISON TO NUMERICAL RESULTS 

Our goal is to understand the effects of nonlinearities on the baryon acoustic peak in the correlation function. We 
have shown that we can derive an expression for the first nonlinear term in configuration space from the Zel'dovich 
approximation. For the real-sp ace correlation function, we a re able to check our analytical expression against the 



Zel'dovich power spectrum from Grinstein & Wise (1987) and Valageas (2011), which is calculated using an entirely 
different technique. This provides analytical validation of our configuration-space approach to perturbation theory. 

However, in the future we would like to extend this calculation to redshift space and to higher orders. For these 
applications, the Fourier equivalents do not exist, so we would have no analytical validation of the results. Therefore, 
we need a numerical method that can be used to test these future results. Here, we present the results from a set of 
Zel'dovich simulations, and show that they validate our configuration-space result, and can therefore be used to test 
future analytical results in redshift space and at higher orders. 

A set of 100 Zel'dovich simulations, written in Python, was run with the cosmological parameters Q,^ — 0.71, 
rim = 0.29, = 0.045, h — 0.7, and <Jg,{z = 0) = 0.89. Each simulation was a 1 Gpc/h box with 512'^ particles, 
and a cell length of 3.9 Mpc/h. We computed the density at 4 different redshifts in each simulation from z — 10 and 
z = using a cloud-in-cell interpolation scheme. The density was then used to compute the correlation function at 
each time step. The average of the correlation functions over the 100 simulations allows us to see the behavior of the 
nonlinear term over time. 

We expect the results from the simulations to agree with the analytical expression at high redshift, although the 
result will be noisy due to the small value of the nonlinear contribution at early times. There will be a growing 
deviation from the analytical result with decreasing redshift coming from the next higher order term, proportional to 
D^. At z = 1, this higher order term will be smaller than 1% of the value at the acoustic peak. At z = 0, however, this 
term will cause a deviation on the order of 5% at the BAO scale, and thus we do not expect the analytical expression 
to approximate well the behavior at z = 0. This is s imilar to the results from 1-loop standard perturbation theory as 
compared to N-body simulations ( Matsubara||2008b ) , where the higher order terms have the greatest contribution at 
low redshift. 

Figure n] shows the averaged correlation functions at redshifts 10, 1, and 0, normalized by for comparison. We 
compare these to the theoretical predictions from the Zel'dovich approximation for the linear and nonlinear correlation 
functions at z = 0. At z = 10, the correlation function is well-approximated by linear theory. At z = 1, the behavior 
of the peak is closer to that of the nonlinear correlation function. At z = 0, higher order terms become non-negligible, 
so the first nonlinear term is not sufficient to fully describe the correlation function. As expected from perturbation 
theory, the acoustic peak is noticeably dampened at later times. 

Figure [2] shows only the nonlinear contribution to the correlation function at z = 10, z = 5, z = 1, and z = 0, 
again normalized by for comparison. We also show the prediction from perturbation theory for the nonlinear term. 
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r [Mpc/h] 

Fig. 1. — The full correlation function. The dashed lines represent the theoretical predictions for the linear (dark green) and nonlinear 
(purple) correlation functions at z=0 from the Zel'dovich approximation. The solid lines are the correlation functions, normalized by D^, 
from 100 realizations. At z = 10, there is good agreement with linear theory. At 2 = 1 the behavior of the peak is well-described by the 
nonlinear correlation function. At 2 = the first nonlinear term is not sufficient to describe the deviation from linear theory, as higher 
order terms become non-negligible. 
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r [Mpc/h] 

Fig. 2. — The nonlinear correlation function. The dashed line represents the nonlinear term from the Zel'dovich approximation. The 
solid lines are the nonlinear term from Zel'dovich simulations, averaged over 100 realizations. At 2 = 10 there is good agreement with the 
analytical prediction. At 2 = 0, higher order terms have a non-negligible effect. 

There is overall agreement with theory at scales larger than about 20 Mpc/h. Below this scale, we are limited by the 
resolution of our simulation. At scales greater than about 140 Mpc/h, we see deviation in the nonlinear term from 
the theoretical prediction. We believe that this is also an artifact of the simulation, related to aliasing of the particle 
window function. The deviation is amplified in this figure due to the multiplication with a very large value of . On 
the scale of the acoustic peak, there is good agreement between the analytical expression and the simulations. At 
z = we again see that the first nonlinear term is not sufficient to accurately describe the peak, because the higher 
order terms become non-negligible. 
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5. CONCLUSION 

We have shown that the first nonhnear term in the correlation function can be directly calculated in configuration 

space from the Zel'dovich approximation, and that this term is in agreement with both the analytical Fourier-space 
result and with numerical simulations. In configuration space, this nonlinear term has a simple form the sum of 
products of a broader class of linear correlation functions, Cn{r). Higher order terms arc in principle straightforward 
to compute, involving terms with products of more than two of these functions. The simple form of these higher order 
terms is one advantage to developing perturbation theory in configuration space. 

Another advantage is that we can straightforwardly extend the real-space correlation function into redshift space. 
This allows us to fully characterize the nonlinear contribution to the acoustic peak in redshift space. In a future paper, 
we will show how this calculation, which for the Zel'dovich approximation has not been possible in Fourier space, can 
be done easily using our configuration-space approach. 
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